TimeSeries
Sergio Solano
9 de marzo de 2017
Pruebas series de tiempo de Time Series Analysis and its applications
# Plot AirPassengers
plot(AirPassengers)
# View the start and end dates of AirPassengers
start(AirPassengers)## [1] 1949 1end(AirPassengers)## [1] 1960 12# Use time(), deltat(), frequency(), and cycle() with AirPassengers
time(AirPassengers)## Jan Feb Mar Apr May Jun Jul
## 1949 1949.000 1949.083 1949.167 1949.250 1949.333 1949.417 1949.500
## 1950 1950.000 1950.083 1950.167 1950.250 1950.333 1950.417 1950.500
## 1951 1951.000 1951.083 1951.167 1951.250 1951.333 1951.417 1951.500
## 1952 1952.000 1952.083 1952.167 1952.250 1952.333 1952.417 1952.500
## 1953 1953.000 1953.083 1953.167 1953.250 1953.333 1953.417 1953.500
## 1954 1954.000 1954.083 1954.167 1954.250 1954.333 1954.417 1954.500
## 1955 1955.000 1955.083 1955.167 1955.250 1955.333 1955.417 1955.500
## 1956 1956.000 1956.083 1956.167 1956.250 1956.333 1956.417 1956.500
## 1957 1957.000 1957.083 1957.167 1957.250 1957.333 1957.417 1957.500
## 1958 1958.000 1958.083 1958.167 1958.250 1958.333 1958.417 1958.500
## 1959 1959.000 1959.083 1959.167 1959.250 1959.333 1959.417 1959.500
## 1960 1960.000 1960.083 1960.167 1960.250 1960.333 1960.417 1960.500
## Aug Sep Oct Nov Dec
## 1949 1949.583 1949.667 1949.750 1949.833 1949.917
## 1950 1950.583 1950.667 1950.750 1950.833 1950.917
## 1951 1951.583 1951.667 1951.750 1951.833 1951.917
## 1952 1952.583 1952.667 1952.750 1952.833 1952.917
## 1953 1953.583 1953.667 1953.750 1953.833 1953.917
## 1954 1954.583 1954.667 1954.750 1954.833 1954.917
## 1955 1955.583 1955.667 1955.750 1955.833 1955.917
## 1956 1956.583 1956.667 1956.750 1956.833 1956.917
## 1957 1957.583 1957.667 1957.750 1957.833 1957.917
## 1958 1958.583 1958.667 1958.750 1958.833 1958.917
## 1959 1959.583 1959.667 1959.750 1959.833 1959.917
## 1960 1960.583 1960.667 1960.750 1960.833 1960.917deltat(AirPassengers)## [1] 0.08333333frequency(AirPassengers)## [1] 12cycle(AirPassengers)## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1949 1 2 3 4 5 6 7 8 9 10 11 12
## 1950 1 2 3 4 5 6 7 8 9 10 11 12
## 1951 1 2 3 4 5 6 7 8 9 10 11 12
## 1952 1 2 3 4 5 6 7 8 9 10 11 12
## 1953 1 2 3 4 5 6 7 8 9 10 11 12
## 1954 1 2 3 4 5 6 7 8 9 10 11 12
## 1955 1 2 3 4 5 6 7 8 9 10 11 12
## 1956 1 2 3 4 5 6 7 8 9 10 11 12
## 1957 1 2 3 4 5 6 7 8 9 10 11 12
## 1958 1 2 3 4 5 6 7 8 9 10 11 12
## 1959 1 2 3 4 5 6 7 8 9 10 11 12
## 1960 1 2 3 4 5 6 7 8 9 10 11 12# Simulate a WN model with list(order = c(0, 0, 0))
white_noise <- arima.sim(model = list(order = c(0,0,0)), n = 50)
# Plot your white_noise data
ts.plot(white_noise)
# Simulate from the WN model with: mean = 100, sd = 10
white_noise_2 <- arima.sim(model = list(order = c(0,0,0)), n = 50, mean = 100, sd = 10)
# Plot your white_noise_2 data
ts.plot(white_noise_2)
Los cambios en una serie de tiempo de Random Walk siguen un comportamiento de White noise.
# Generate a RW model using arima.sim
random_walk <- arima.sim(model = list(order = c(0, 1, 0)) , n = 100)
# Plot random_walk
ts.plot (random_walk)
# Calculate the first difference series
random_walk_diff <- diff(random_walk)
# Plot random_walk_diff
ts.plot(random_walk_diff)
# RANDOM WALK WITH DRIFT
# Generate a RW model with a drift uing arima.sim
rw_drift <- arima.sim(model = list(order = c(0, 1, 0)), n = 100, mean = 1)
# Plot rw_drift
ts.plot(rw_drift)
# Calculate the first difference series
rw_drift_diff <- diff(rw_drift)
# Plot rw_drift_diff
ts.plot(rw_drift_diff)
random_walk <- rw_drift
# Difference your random_walk data
rw_diff <- diff(random_walk)
# Plot rw_diff
plot.ts(rw_diff)
# Now fit the WN model to the differenced data
model_wn <- arima(rw_diff,order = c(0, 0, 0))
# Store the value of the estimated time trend (intercept)
int_wn <- model_wn$coef
# Plot the original random_walk data
ts.plot(random_walk)
# Use abline(0, ...) to add time trend to the figure
abline(0,int_wn)
# RANDOM WALK CON Y SIN DRIFT (Estationarity)
#
# The white noise (WN) and random walk (RW) models are very closely related. However, only the RW is always non-stationary, both with and without a drift term. This is a simulation exercise to highlight the differences.
#
# Recall that if we start with a mean zero WN process and compute its running or cumulative sum, the result is a RW process. The cumsum() function will make this transformation for you. Similarly, if we create a WN process, but change its mean from zero, and then compute its cumulative sum, the result is a RW process with a drift.
# Use arima.sim() to generate WN data
white_noise <- arima.sim(model = list(order = c(0, 0, 0)), n = 100)
# Use cumsum() to convert your WN data to RW
random_walk <- cumsum(white_noise)
# Use arima.sim() to generate WN drift data
wn_drift <- arima.sim(model = list(order = c(0, 0, 0)),mean=0.4, n = 100)
# Use cumsum() to convert your WN drift data to RW
rw_drift <- cumsum(wn_drift)
# Plot all four data objects
plot.ts(cbind(white_noise, random_walk, wn_drift, rw_drift))
#
# # Generate means from eu_percentreturns
# colMeans(eu_percentreturns)
#
# # Use apply to calculate sample variance from eu_percentreturns
# apply(eu_percentreturns, MARGIN = 2, FUN = var)
#
# # Use apply to calculate standard deviation from eu_percentreturns
# apply(eu_percentreturns, MARGIN = 2, FUN = sd)
#
# # Display a histogram of percent returns for each index
# par(mfrow = c(2,2))
# apply(eu_percentreturns, MARGIN = 2, FUN = hist, main = "", xlab = "Percentage Return")
#
# # Display normal quantile plots of percent returns for each index
# par(mfrow = c(2,2))
# apply(eu_percentreturns, MARGIN = 2, FUN = qqnorm, main = "")
# qqline(eu_percentreturns)
# pairs(eu_stocks)
# MODELOS AUTOREGRESIVOS:
# Simulate an AR model with 0.5 slope
x <- arima.sim(model = list(ar = 0.5), n = 100)
# Simulate an AR model with 0.9 slope
y <- arima.sim(model = list(ar = 0.9), n = 100)
# Simulate an AR model with -0.75 slope
z <- arima.sim(model = list(ar = -0.75), n = 100)
# Plot your simulated data
plot.ts(cbind(x, y, z))
###### COMPARAR MODELOS AR con RandomWalks
# Simulate and plot AR model with slope 0.9
x <- arima.sim(model = list(ar = 0.9), n = 200)
ts.plot(x)
acf(x)
# Simulate and plot AR model with slope 0.98
y <- arima.sim(model = list(ar = 0.98), n = 200)
ts.plot(y)
acf(y)
# Simulate and plot RW model
z <- arima.sim(model = list(order = c(0, 1, 0)), n = 200)
ts.plot(z)
acf(z)
# # AJUSTAR UN AR
# # Fit an AR model to Nile
# AR_fit <-arima(Nile, order = c(1,0,0))
# print(AR_fit)
#
# # Use predict() to make a 1-step forecast
# predict_AR <- predict(AR_fit)
#
# # Obtain the 1-step forecast using $pred[1]
# predict_AR$pred[1]
#
# # Use predict to make 1-step through 10-step forecasts
# predict(AR_fit, n.ahead = 10)
#
# # Run to plot the Nile series plus the forecast and 95% prediction intervals
# ts.plot(Nile, xlim = c(1871, 1980))
# AR_forecast <- predict(AR_fit, n.ahead = 10)$pred
# AR_forecast_se <- predict(AR_fit, n.ahead = 10)$se
# points(AR_forecast, type = "l", col = 2)
# points(AR_forecast - 2*AR_forecast_se, type = "l", col = 2, lty = 2)
# points(AR_forecast + 2*AR_forecast_se, type = "l", col = 2, lty = 2)
# # AJUSTAR UN MA
# Generate MA model with slope 0.5
x <- arima.sim(model = list(ma = 0.5), n = 100)
# Generate MA model with slope 0.9
y <- x <- arima.sim(model = list(ma = 0.9), n = 100)
# Generate MA model with slope -0.5
z <- x <- arima.sim(model = list(ma = -0.5), n = 100)
# Plot all three models together
plot.ts(cbind(x, y, z))
# # # Ajustando un MA
#
# # Make a 1-step forecast based on MA
# predict_MA <- predict(MA, n.ahead = 1)
#
# # Obtain the 1-step forecast using $pred[1]
# predict_MA$pred[1]
#
# # Make a 1-step through 10-step forecast based on MA
# predict(MA, n.ahead = 10)
#
# # Plot the Nile series plus the forecast and 95% prediction intervals
# ts.plot(Nile, xlim = c(1871, 1980))
# MA_forecasts <- predict(MA, n.ahead = 10)$pred
# MA_forecast_se <- predict(MA, n.ahead = 10)$se
# points(MA_forecasts, type = "l", col = 2)
# points(MA_forecasts - 2*MA_forecast_se, type = "l", col = 2, lty = 2)
# points(MA_forecasts + 2*MA_forecast_se, type = "l", col = 2, lty = 2)Pruebas series de tiempo de Time Series Analysis and its applications course
# # R Arima base
# View a detailed description of AirPassengers
help(AirPassengers)
# Plot AirPassengers
ts.plot(AirPassengers)
# Plot the DJIA daily closings
ts.plot(djia$Close)
# Plot the Southern Oscillation Index
plot(soi)
# Plot globtemp and detrended globtemp
par(mfrow = c(2,1))
plot(globtemp)
plot(diff(globtemp))
# Plot cmort and detrended cmort
par(mfrow = c(2,1))
plot(cmort)
plot(diff(cmort))
# Plot GNP series (gnp) and its growth rate
par(mfrow = c(2,1))
plot(gnp)
plot(diff(log(gnp)))
# Plot DJIA closings (djia$Close) and its returns
par(mfrow = c(2,1))
plot(djia$Close)
plot(diff(log(djia$Close)))
##SIMULAR MA & AR
# Generate and plot white noise
WN <- arima.sim(model=list(order=c(0,0,0)),n=100)
plot(WN)
# Generate and plot an MA(1) with parameter .9
MA <- arima.sim(model=list(order=c(0,0,1), ma = 0.9),n=100)
plot(MA)
# Generate and plot an AR(2) with parameters 1.5 and -.75
AR <- arima.sim(model=list(order=c(2,0,0), ar =c(1.5,-0.75)),n=100)
plot(AR)
# # SIMULAR Y AJUSTAR AR
# Generate 100 observations from the AR(1) model
x <- arima.sim(model = list(order = c(1, 0, 0), ar = .9), n = 100)
# Plot the generated data
plot(x)
# Plot the sample P/ACF pair
acf2(x)
## ACF PACF
## [1,] 0.92 0.92
## [2,] 0.86 0.05
## [3,] 0.79 -0.04
## [4,] 0.74 0.06
## [5,] 0.69 0.01
## [6,] 0.65 0.03
## [7,] 0.62 0.06
## [8,] 0.56 -0.18
## [9,] 0.53 0.14
## [10,] 0.49 -0.10
## [11,] 0.46 0.06
## [12,] 0.41 -0.13
## [13,] 0.36 -0.09
## [14,] 0.31 0.00
## [15,] 0.26 -0.03
## [16,] 0.23 0.04
## [17,] 0.19 -0.03
## [18,] 0.17 0.03
## [19,] 0.15 0.05
## [20,] 0.13 -0.03# Fit an AR(1) to the data and examine the -table
sarima(x, 1, 0, 0)## initial value 0.959788
## iter 2 value 0.007398
## iter 3 value 0.007395
## iter 4 value 0.007392
## iter 5 value 0.007373
## iter 6 value 0.007366
## iter 7 value 0.007364
## iter 8 value 0.007364
## iter 9 value 0.007364
## iter 10 value 0.007364
## iter 11 value 0.007364
## iter 12 value 0.007364
## iter 13 value 0.007363
## iter 14 value 0.007363
## iter 15 value 0.007363
## iter 16 value 0.007363
## iter 17 value 0.007363
## iter 17 value 0.007363
## iter 17 value 0.007363
## final value 0.007363
## converged
## initial value 0.012435
## iter 2 value 0.012362
## iter 3 value 0.012164
## iter 4 value 0.012152
## iter 5 value 0.012141
## iter 6 value 0.012114
## iter 7 value 0.012113
## iter 8 value 0.012108
## iter 9 value 0.012108
## iter 10 value 0.012106
## iter 11 value 0.012106
## iter 11 value 0.012106
## final value 0.012106
## converged
## $fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), xreg = xmean, include.mean = FALSE, optim.control = list(trace = trc,
## REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 xmean
## 0.9147 -0.9672
## s.e. 0.0367 1.0668
##
## sigma^2 estimated as 1.006: log likelihood = -143.1, aic = 292.21
##
## $degrees_of_freedom
## [1] 98
##
## $ttable
## Estimate SE t.value p.value
## ar1 0.9147 0.0367 24.8895 0.0000
## xmean -0.9672 1.0668 -0.9066 0.3668
##
## $AIC
## [1] 1.046097
##
## $AICc
## [1] 1.068597
##
## $BIC
## [1] 0.09820056# # aJUSTAR MA(2)
# astsa is preloaded
x <- arima.sim(model = list(order = c(2, 0, 0), ar = c(1.5, -.75)), n = 200)
# Plot x
plot(x)
# Plot the sample P/ACF of x
acf2(x)
## ACF PACF
## [1,] 0.84 0.84
## [2,] 0.49 -0.71
## [3,] 0.09 -0.07
## [4,] -0.23 -0.04
## [5,] -0.43 -0.03
## [6,] -0.47 -0.06
## [7,] -0.41 -0.03
## [8,] -0.25 0.12
## [9,] -0.03 0.12
## [10,] 0.18 0.05
## [11,] 0.34 -0.05
## [12,] 0.39 0.01
## [13,] 0.32 -0.01
## [14,] 0.17 -0.09
## [15,] -0.03 -0.08
## [16,] -0.22 -0.05
## [17,] -0.36 -0.06
## [18,] -0.40 0.04
## [19,] -0.34 -0.01
## [20,] -0.21 -0.07
## [21,] -0.05 -0.07
## [22,] 0.07 -0.09
## [23,] 0.15 0.07
## [24,] 0.17 -0.06
## [25,] 0.13 -0.06# Fit an AR(2) to the data and examine the t-table
sarima(x,2,0,0)## initial value 1.003984
## iter 2 value 0.849174
## iter 3 value 0.465756
## iter 4 value 0.241604
## iter 5 value 0.103569
## iter 6 value 0.005471
## iter 7 value -0.013685
## iter 8 value -0.019326
## iter 9 value -0.019331
## iter 10 value -0.019348
## iter 11 value -0.019357
## iter 12 value -0.019359
## iter 13 value -0.019359
## iter 13 value -0.019359
## iter 13 value -0.019359
## final value -0.019359
## converged
## initial value -0.008804
## iter 2 value -0.008826
## iter 3 value -0.008840
## iter 4 value -0.008841
## iter 5 value -0.008841
## iter 6 value -0.008841
## iter 7 value -0.008841
## iter 7 value -0.008841
## iter 7 value -0.008841
## final value -0.008841
## converged
## $fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), xreg = xmean, include.mean = FALSE, optim.control = list(trace = trc,
## REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 ar2 xmean
## 1.4695 -0.7449 0.2919
## s.e. 0.0472 0.0471 0.2529
##
## sigma^2 estimated as 0.9686: log likelihood = -282.02, aic = 572.04
##
## $degrees_of_freedom
## [1] 197
##
## $ttable
## Estimate SE t.value p.value
## ar1 1.4695 0.0472 31.1444 0.0000
## ar2 -0.7449 0.0471 -15.7996 0.0000
## xmean 0.2919 0.2529 1.1543 0.2498
##
## $AIC
## [1] 0.9980474
##
## $AICc
## [1] 1.009073
##
## $BIC
## [1] 0.04752219# # # SIMULANDO ARIMAS
# Fit an ARIMA(0,1,2) to globtemp and check the fit
sarima(globtemp, 0, 1, 2)## initial value -2.220513
## iter 2 value -2.294887
## iter 3 value -2.307682
## iter 4 value -2.309170
## iter 5 value -2.310360
## iter 6 value -2.311251
## iter 7 value -2.311636
## iter 8 value -2.311648
## iter 9 value -2.311649
## iter 9 value -2.311649
## iter 9 value -2.311649
## final value -2.311649
## converged
## initial value -2.310187
## iter 2 value -2.310197
## iter 3 value -2.310199
## iter 4 value -2.310201
## iter 5 value -2.310202
## iter 5 value -2.310202
## iter 5 value -2.310202
## final value -2.310202
## converged
## $fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), xreg = constant, optim.control = list(trace = trc, REPORT = 1,
## reltol = tol))
##
## Coefficients:
## ma1 ma2 constant
## -0.3984 -0.2173 0.0072
## s.e. 0.0808 0.0768 0.0033
##
## sigma^2 estimated as 0.00982: log likelihood = 120.32, aic = -232.64
##
## $degrees_of_freedom
## [1] 133
##
## $ttable
## Estimate SE t.value p.value
## ma1 -0.3984 0.0808 -4.9313 0.0000
## ma2 -0.2173 0.0768 -2.8303 0.0054
## constant 0.0072 0.0033 2.1463 0.0337
##
## $AIC
## [1] -3.579224
##
## $AICc
## [1] -3.562273
##
## $BIC
## [1] -4.514974# Forecast data 35 years into the future
sarima.for (globtemp,n.ahead = 35, p = 0, d = 1, q = 2) ## $pred
## Time Series:
## Start = 2016
## End = 2050
## Frequency = 1
## [1] 0.7995567 0.7745381 0.7816919 0.7888457 0.7959996 0.8031534 0.8103072
## [8] 0.8174611 0.8246149 0.8317688 0.8389226 0.8460764 0.8532303 0.8603841
## [15] 0.8675379 0.8746918 0.8818456 0.8889995 0.8961533 0.9033071 0.9104610
## [22] 0.9176148 0.9247687 0.9319225 0.9390763 0.9462302 0.9533840 0.9605378
## [29] 0.9676917 0.9748455 0.9819994 0.9891532 0.9963070 1.0034609 1.0106147
##
## $se
## Time Series:
## Start = 2016
## End = 2050
## Frequency = 1
## [1] 0.09909556 0.11564576 0.12175580 0.12757353 0.13313729 0.13847769
## [7] 0.14361964 0.14858376 0.15338730 0.15804492 0.16256915 0.16697084
## [13] 0.17125943 0.17544322 0.17952954 0.18352490 0.18743511 0.19126540
## [19] 0.19502047 0.19870459 0.20232164 0.20587515 0.20936836 0.21280424
## [25] 0.21618551 0.21951471 0.22279416 0.22602604 0.22921235 0.23235497
## [31] 0.23545565 0.23851603 0.24153763 0.24452190 0.24747019### SERIES DE TIEMPO ESTACIONALES
# Plot unemp
plot(unemp)
# Difference your data and plot it
d_unemp <- diff(unemp)
plot(d_unemp)
# Seasonally difference d_unemp and plot it
dd_unemp <- diff(d_unemp, lag = 12)
plot(dd_unemp)
# Plot P/ACF pair of fully differenced data to lag 60
dd_unemp <- diff(diff(unemp), lag = 12)
acf2(dd_unemp, max.lag = 60)
## ACF PACF
## [1,] 0.21 0.21
## [2,] 0.33 0.29
## [3,] 0.15 0.05
## [4,] 0.17 0.05
## [5,] 0.10 0.01
## [6,] 0.06 -0.02
## [7,] -0.06 -0.12
## [8,] -0.02 -0.03
## [9,] -0.09 -0.05
## [10,] -0.17 -0.15
## [11,] -0.08 0.02
## [12,] -0.48 -0.43
## [13,] -0.18 -0.02
## [14,] -0.16 0.15
## [15,] -0.11 0.03
## [16,] -0.15 -0.04
## [17,] -0.09 -0.01
## [18,] -0.09 0.00
## [19,] 0.03 0.01
## [20,] -0.01 0.01
## [21,] 0.02 -0.01
## [22,] -0.02 -0.16
## [23,] 0.01 0.01
## [24,] -0.02 -0.27
## [25,] 0.09 0.05
## [26,] -0.05 -0.01
## [27,] -0.01 -0.05
## [28,] 0.03 0.05
## [29,] 0.08 0.09
## [30,] 0.01 -0.04
## [31,] 0.03 0.02
## [32,] -0.05 -0.07
## [33,] 0.01 -0.01
## [34,] 0.02 -0.08
## [35,] -0.06 -0.08
## [36,] -0.02 -0.23
## [37,] -0.12 -0.08
## [38,] 0.01 0.06
## [39,] -0.03 -0.07
## [40,] -0.03 -0.01
## [41,] -0.10 0.03
## [42,] -0.02 -0.03
## [43,] -0.13 -0.11
## [44,] 0.00 -0.04
## [45,] -0.06 0.01
## [46,] 0.01 0.00
## [47,] 0.02 -0.03
## [48,] 0.11 -0.04
## [49,] 0.13 0.02
## [50,] 0.10 0.03
## [51,] 0.07 -0.05
## [52,] 0.10 0.02
## [53,] 0.12 0.02
## [54,] 0.06 -0.08
## [55,] 0.14 0.00
## [56,] 0.05 -0.03
## [57,] 0.04 -0.07
## [58,] 0.04 0.05
## [59,] 0.07 0.04
## [60,] -0.03 -0.04# Fit an appropriate model
sarima(unemp,2,1,0,0,1,1, S=12)## initial value 3.340809
## iter 2 value 3.105512
## iter 3 value 3.086631
## iter 4 value 3.079778
## iter 5 value 3.069447
## iter 6 value 3.067659
## iter 7 value 3.067426
## iter 8 value 3.067418
## iter 8 value 3.067418
## final value 3.067418
## converged
## initial value 3.065481
## iter 2 value 3.065478
## iter 3 value 3.065477
## iter 3 value 3.065477
## iter 3 value 3.065477
## final value 3.065477
## converged
## $fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), include.mean = !no.constant, optim.control = list(trace = trc,
## REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 ar2 sma1
## 0.1351 0.2464 -0.6953
## s.e. 0.0513 0.0515 0.0381
##
## sigma^2 estimated as 449.6: log likelihood = -1609.91, aic = 3227.81
##
## $degrees_of_freedom
## [1] 369
##
## $ttable
## Estimate SE t.value p.value
## ar1 0.1351 0.0513 2.6326 0.0088
## ar2 0.2464 0.0515 4.7795 0.0000
## sma1 -0.6953 0.0381 -18.2362 0.0000
##
## $AIC
## [1] 7.12457
##
## $AICc
## [1] 7.130239
##
## $BIC
## [1] 6.156174# Plot differenced chicken
plot(diff(chicken))
# Plot P/ACF pair of differenced data to lag 60
acf2(diff(chicken), max.lag = 60)
## ACF PACF
## [1,] 0.72 0.72
## [2,] 0.39 -0.29
## [3,] 0.09 -0.14
## [4,] -0.07 0.03
## [5,] -0.16 -0.10
## [6,] -0.20 -0.06
## [7,] -0.27 -0.19
## [8,] -0.23 0.12
## [9,] -0.11 0.10
## [10,] 0.09 0.16
## [11,] 0.26 0.09
## [12,] 0.33 0.00
## [13,] 0.20 -0.22
## [14,] 0.07 0.03
## [15,] -0.03 0.03
## [16,] -0.10 -0.11
## [17,] -0.19 -0.09
## [18,] -0.25 0.01
## [19,] -0.29 -0.03
## [20,] -0.20 0.07
## [21,] -0.08 -0.04
## [22,] 0.08 0.06
## [23,] 0.16 -0.05
## [24,] 0.18 0.02
## [25,] 0.08 -0.14
## [26,] -0.06 -0.19
## [27,] -0.21 -0.13
## [28,] -0.31 -0.06
## [29,] -0.40 -0.08
## [30,] -0.40 -0.05
## [31,] -0.33 0.01
## [32,] -0.18 0.03
## [33,] 0.02 0.10
## [34,] 0.20 0.02
## [35,] 0.30 -0.01
## [36,] 0.35 0.09
## [37,] 0.26 -0.12
## [38,] 0.13 0.01
## [39,] -0.02 -0.01
## [40,] -0.14 -0.05
## [41,] -0.23 0.02
## [42,] -0.21 0.12
## [43,] -0.18 -0.05
## [44,] -0.11 -0.13
## [45,] -0.03 -0.07
## [46,] 0.08 0.01
## [47,] 0.21 0.14
## [48,] 0.33 0.05
## [49,] 0.26 -0.20
## [50,] 0.12 -0.01
## [51,] -0.01 0.07
## [52,] -0.11 -0.04
## [53,] -0.13 0.02
## [54,] -0.09 0.00
## [55,] -0.09 -0.08
## [56,] -0.06 0.03
## [57,] 0.03 0.04
## [58,] 0.17 0.00
## [59,] 0.29 0.01
## [60,] 0.32 0.03# Fit ARIMA(2,1,0) to chicken - not so good
sarima(chicken, 2,1,0)## initial value 0.001863
## iter 2 value -0.156034
## iter 3 value -0.359181
## iter 4 value -0.424164
## iter 5 value -0.430212
## iter 6 value -0.432744
## iter 7 value -0.432747
## iter 8 value -0.432749
## iter 9 value -0.432749
## iter 10 value -0.432751
## iter 11 value -0.432752
## iter 12 value -0.432752
## iter 13 value -0.432752
## iter 13 value -0.432752
## iter 13 value -0.432752
## final value -0.432752
## converged
## initial value -0.420883
## iter 2 value -0.420934
## iter 3 value -0.420936
## iter 4 value -0.420937
## iter 5 value -0.420937
## iter 6 value -0.420937
## iter 6 value -0.420937
## iter 6 value -0.420937
## final value -0.420937
## converged
## $fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), xreg = constant, optim.control = list(trace = trc, REPORT = 1,
## reltol = tol))
##
## Coefficients:
## ar1 ar2 constant
## 0.9494 -0.3069 0.2632
## s.e. 0.0717 0.0718 0.1362
##
## sigma^2 estimated as 0.4286: log likelihood = -178.64, aic = 365.28
##
## $degrees_of_freedom
## [1] 177
##
## $ttable
## Estimate SE t.value p.value
## ar1 0.9494 0.0717 13.2339 0.0000
## ar2 -0.3069 0.0718 -4.2723 0.0000
## constant 0.2632 0.1362 1.9328 0.0549
##
## $AIC
## [1] 0.1861622
##
## $AICc
## [1] 0.1985432
##
## $BIC
## [1] -0.7606218# Fit SARIMA(2,1,0,1,0,0,12) to chicken - that works
sarima(chicken,2,1,0,1,0,0,12)## initial value 0.015039
## iter 2 value -0.226398
## iter 3 value -0.412955
## iter 4 value -0.460882
## iter 5 value -0.470787
## iter 6 value -0.471082
## iter 7 value -0.471088
## iter 8 value -0.471090
## iter 9 value -0.471092
## iter 10 value -0.471095
## iter 11 value -0.471095
## iter 12 value -0.471096
## iter 13 value -0.471096
## iter 14 value -0.471096
## iter 15 value -0.471097
## iter 16 value -0.471097
## iter 16 value -0.471097
## iter 16 value -0.471097
## final value -0.471097
## converged
## initial value -0.473585
## iter 2 value -0.473664
## iter 3 value -0.473721
## iter 4 value -0.473823
## iter 5 value -0.473871
## iter 6 value -0.473885
## iter 7 value -0.473886
## iter 8 value -0.473886
## iter 8 value -0.473886
## iter 8 value -0.473886
## final value -0.473886
## converged
## $fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), xreg = constant, optim.control = list(trace = trc, REPORT = 1,
## reltol = tol))
##
## Coefficients:
## ar1 ar2 sar1 constant
## 0.9154 -0.2494 0.3237 0.2353
## s.e. 0.0733 0.0739 0.0715 0.1973
##
## sigma^2 estimated as 0.3828: log likelihood = -169.16, aic = 348.33
##
## $degrees_of_freedom
## [1] 176
##
## $ttable
## Estimate SE t.value p.value
## ar1 0.9154 0.0733 12.4955 0.0000
## ar2 -0.2494 0.0739 -3.3728 0.0009
## sar1 0.3237 0.0715 4.5238 0.0000
## constant 0.2353 0.1973 1.1923 0.2347
##
## $AIC
## [1] 0.0842377
##
## $AICc
## [1] 0.09726452
##
## $BIC
## [1] -0.8448077###### Natalidad US
# Plot P/ACF to lag 60 of differenced data
d_birth <- diff(birth)
acf2(d_birth)
## ACF PACF
## [1,] -0.32 -0.32
## [2,] 0.16 0.06
## [3,] -0.08 -0.01
## [4,] -0.19 -0.25
## [5,] 0.09 -0.03
## [6,] -0.28 -0.26
## [7,] 0.06 -0.17
## [8,] -0.19 -0.29
## [9,] -0.05 -0.35
## [10,] 0.17 -0.16
## [11,] -0.26 -0.59
## [12,] 0.82 0.57
## [13,] -0.28 0.13
## [14,] 0.17 0.11
## [15,] -0.07 0.13
## [16,] -0.18 0.09
## [17,] 0.08 0.00
## [18,] -0.28 0.00
## [19,] 0.07 0.05
## [20,] -0.18 0.04
## [21,] -0.05 -0.07
## [22,] 0.16 -0.10
## [23,] -0.24 -0.20
## [24,] 0.78 0.19
## [25,] -0.27 0.01
## [26,] 0.19 0.05
## [27,] -0.08 0.07
## [28,] -0.17 0.07
## [29,] 0.07 -0.02
## [30,] -0.29 -0.06# Plot P/ACF to lag 60 of seasonal differenced data
dd_birth <- diff(d_birth, lag = 12)
acf2(dd_birth)
## ACF PACF
## [1,] -0.30 -0.30
## [2,] -0.09 -0.20
## [3,] -0.09 -0.21
## [4,] 0.00 -0.14
## [5,] 0.07 -0.03
## [6,] 0.03 0.02
## [7,] -0.07 -0.06
## [8,] -0.04 -0.08
## [9,] 0.11 0.06
## [10,] 0.04 0.08
## [11,] 0.13 0.23
## [12,] -0.43 -0.32
## [13,] 0.14 -0.06
## [14,] -0.01 -0.13
## [15,] 0.03 -0.13
## [16,] 0.01 -0.11
## [17,] 0.02 0.02
## [18,] 0.00 0.06
## [19,] 0.03 0.04
## [20,] -0.07 -0.10
## [21,] -0.01 0.02
## [22,] 0.00 0.00
## [23,] 0.06 0.17
## [24,] -0.01 -0.13
## [25,] -0.12 -0.14
## [26,] 0.17 0.07
## [27,] -0.04 -0.04
## [28,] 0.03 -0.02
## [29,] -0.05 0.02# Fit SARIMA(0,1,1)x(0,1,1)_12. What happens?
sarima(birth,0,1,1,0,1,1,12)## initial value 2.219164
## iter 2 value 2.013310
## iter 3 value 1.988107
## iter 4 value 1.980026
## iter 5 value 1.967594
## iter 6 value 1.965384
## iter 7 value 1.965049
## iter 8 value 1.964993
## iter 9 value 1.964992
## iter 9 value 1.964992
## iter 9 value 1.964992
## final value 1.964992
## converged
## initial value 1.951264
## iter 2 value 1.945867
## iter 3 value 1.945729
## iter 4 value 1.945723
## iter 5 value 1.945723
## iter 5 value 1.945723
## iter 5 value 1.945723
## final value 1.945723
## converged
## $fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), include.mean = !no.constant, optim.control = list(trace = trc,
## REPORT = 1, reltol = tol))
##
## Coefficients:
## ma1 sma1
## -0.4734 -0.7861
## s.e. 0.0598 0.0451
##
## sigma^2 estimated as 47.4: log likelihood = -1211.28, aic = 2428.56
##
## $degrees_of_freedom
## [1] 371
##
## $ttable
## Estimate SE t.value p.value
## ma1 -0.4734 0.0598 -7.9097 0
## sma1 -0.7861 0.0451 -17.4227 0
##
## $AIC
## [1] 4.869388
##
## $AICc
## [1] 4.874924
##
## $BIC
## [1] 3.890415# Add AR term and conclude
sarima(birth,1,1,1,0,1,1,12)## initial value 2.218186
## iter 2 value 2.032584
## iter 3 value 1.982464
## iter 4 value 1.975643
## iter 5 value 1.971721
## iter 6 value 1.967284
## iter 7 value 1.963840
## iter 8 value 1.961106
## iter 9 value 1.960849
## iter 10 value 1.960692
## iter 11 value 1.960683
## iter 12 value 1.960675
## iter 13 value 1.960672
## iter 13 value 1.960672
## iter 13 value 1.960672
## final value 1.960672
## converged
## initial value 1.940459
## iter 2 value 1.934425
## iter 3 value 1.932752
## iter 4 value 1.931750
## iter 5 value 1.931074
## iter 6 value 1.930882
## iter 7 value 1.930860
## iter 8 value 1.930859
## iter 8 value 1.930859
## final value 1.930859
## converged
## $fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), include.mean = !no.constant, optim.control = list(trace = trc,
## REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 ma1 sma1
## 0.3038 -0.7006 -0.8000
## s.e. 0.0865 0.0604 0.0441
##
## sigma^2 estimated as 45.91: log likelihood = -1205.93, aic = 2419.85
##
## $degrees_of_freedom
## [1] 370
##
## $ttable
## Estimate SE t.value p.value
## ar1 0.3038 0.0865 3.5104 5e-04
## ma1 -0.7006 0.0604 -11.5984 0e+00
## sma1 -0.8000 0.0441 -18.1302 0e+00
##
## $AIC
## [1] 4.842869
##
## $AICc
## [1] 4.848523
##
## $BIC
## [1] 3.87441## FORECASTING SEASONAL DATA
# Fit your previous model to unemp and check the diagnostics
sarima(unemp, 2,1,0,0,1,1,12)## initial value 3.340809
## iter 2 value 3.105512
## iter 3 value 3.086631
## iter 4 value 3.079778
## iter 5 value 3.069447
## iter 6 value 3.067659
## iter 7 value 3.067426
## iter 8 value 3.067418
## iter 8 value 3.067418
## final value 3.067418
## converged
## initial value 3.065481
## iter 2 value 3.065478
## iter 3 value 3.065477
## iter 3 value 3.065477
## iter 3 value 3.065477
## final value 3.065477
## converged
## $fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), include.mean = !no.constant, optim.control = list(trace = trc,
## REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 ar2 sma1
## 0.1351 0.2464 -0.6953
## s.e. 0.0513 0.0515 0.0381
##
## sigma^2 estimated as 449.6: log likelihood = -1609.91, aic = 3227.81
##
## $degrees_of_freedom
## [1] 369
##
## $ttable
## Estimate SE t.value p.value
## ar1 0.1351 0.0513 2.6326 0.0088
## ar2 0.2464 0.0515 4.7795 0.0000
## sma1 -0.6953 0.0381 -18.2362 0.0000
##
## $AIC
## [1] 7.12457
##
## $AICc
## [1] 7.130239
##
## $BIC
## [1] 6.156174# Forecast the data 3 years into the future
sarima.for(unemp,n.ahead = 36,2,1,0,0,1,1,12)## $pred
## Jan Feb Mar Apr May Jun Jul
## 1979 676.4664 685.1172 653.2388 585.6939 553.8813 664.4072 647.0657
## 1980 683.3045 687.7649 654.8658 586.1507 553.9285 664.1108 646.6220
## 1981 682.6406 687.0977 654.1968 585.4806 553.2579 663.4398 645.9508
## Aug Sep Oct Nov Dec
## 1979 611.0828 594.6414 569.3997 587.5801 581.1833
## 1980 610.5345 594.0427 568.7684 586.9320 580.5249
## 1981 609.8632 593.3713 568.0970 586.2606 579.8535
##
## $se
## Jan Feb Mar Apr May Jun Jul
## 1979 21.20465 32.07710 43.70167 53.66329 62.85364 71.12881 78.73590
## 1980 116.99599 124.17344 131.51281 138.60466 145.49706 152.12863 158.52302
## 1981 194.25167 201.10648 208.17066 215.11503 221.96039 228.64285 235.16874
## Aug Sep Oct Nov Dec
## 1979 85.75096 92.28663 98.41329 104.19488 109.67935
## 1980 164.68623 170.63839 176.39520 181.97333 187.38718
## 1981 241.53258 247.74268 253.80549 259.72970 265.52323# Fit the chicken model again and check diagnostics
sarima(chicken,2,1,0,1,0,0,12)
## initial value 0.015039
## iter 2 value -0.226398
## iter 3 value -0.412955
## iter 4 value -0.460882
## iter 5 value -0.470787
## iter 6 value -0.471082
## iter 7 value -0.471088
## iter 8 value -0.471090
## iter 9 value -0.471092
## iter 10 value -0.471095
## iter 11 value -0.471095
## iter 12 value -0.471096
## iter 13 value -0.471096
## iter 14 value -0.471096
## iter 15 value -0.471097
## iter 16 value -0.471097
## iter 16 value -0.471097
## iter 16 value -0.471097
## final value -0.471097
## converged
## initial value -0.473585
## iter 2 value -0.473664
## iter 3 value -0.473721
## iter 4 value -0.473823
## iter 5 value -0.473871
## iter 6 value -0.473885
## iter 7 value -0.473886
## iter 8 value -0.473886
## iter 8 value -0.473886
## iter 8 value -0.473886
## final value -0.473886
## converged
## $fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), xreg = constant, optim.control = list(trace = trc, REPORT = 1,
## reltol = tol))
##
## Coefficients:
## ar1 ar2 sar1 constant
## 0.9154 -0.2494 0.3237 0.2353
## s.e. 0.0733 0.0739 0.0715 0.1973
##
## sigma^2 estimated as 0.3828: log likelihood = -169.16, aic = 348.33
##
## $degrees_of_freedom
## [1] 176
##
## $ttable
## Estimate SE t.value p.value
## ar1 0.9154 0.0733 12.4955 0.0000
## ar2 -0.2494 0.0739 -3.3728 0.0009
## sar1 0.3237 0.0715 4.5238 0.0000
## constant 0.2353 0.1973 1.1923 0.2347
##
## $AIC
## [1] 0.0842377
##
## $AICc
## [1] 0.09726452
##
## $BIC
## [1] -0.8448077# Forecast the chicken data 5 years into the future
sarima.for(chicken,n.ahead = 60,2,1,0,1,0,0,12)## $pred
## Jan Feb Mar Apr May Jun Jul
## 2016
## 2017 110.5358 110.5612 110.5480 110.7055 111.0047 111.1189 111.1552
## 2018 111.8108 111.9782 112.1330 112.3431 112.5991 112.7952 112.9661
## 2019 114.1331 114.3464 114.5556 114.7827 115.0247 115.2473 115.4617
## 2020 116.7942 117.0224 117.2492 117.4819 117.7193 117.9505 118.1790
## 2021 119.5651 119.7980 120.0306 120.2650 120.5010 120.7350 120.9681
## Aug Sep Oct Nov Dec
## 2016 111.0907 110.8740 110.6853 110.5045 110.5527
## 2017 111.1948 111.2838 111.3819 111.4825 111.6572
## 2018 113.1380 113.3260 113.5168 113.7085 113.9242
## 2019 115.6765 115.8965 116.1174 116.3386 116.5675
## 2020 118.4077 118.6380 118.8686 119.0993 119.3326
## 2021
##
## $se
## Jan Feb Mar Apr May Jun
## 2016
## 2017 3.7414959 4.1793190 4.5747009 4.9373266 5.2742129 5.5903499
## 2018 8.2010253 8.5605811 8.9054714 9.2372195 9.5572539 9.8667955
## 2019 12.0038164 12.2921541 12.5738417 12.8492868 13.1188976 13.3830477
## 2020 15.1557253 15.3959082 15.6323906 15.8653300 16.0948844 16.3212022
## 2021 17.8397890 18.0473081 18.2524651 18.4553364 18.6559977 18.8545213
## Jul Aug Sep Oct Nov Dec
## 2016 0.6187194 1.3368594 2.0462419 2.6867986 3.2486625
## 2017 5.8893133 6.2367345 6.6253573 7.0309771 7.4344077 7.8255932
## 2018 10.1668604 10.4736807 10.7857727 11.0980056 11.4063211 11.7085266
## 2019 13.6420693 13.9002670 14.1573839 14.4122197 14.6638269 14.9117124
## 2020 16.5444204 16.7657634 16.9852163 17.2025022 17.4174076 17.6298379
## 2021 19.0509752
Pruebas series de tiempo de Time Series Analysis and its applications text
Página 13 Moving Average
w = rnorm (500,0,1)
v = filter(w, sides = 2, rep(1/3,3))
par(mfrow=c(2,1))
plot.ts(w,main = "white noise")
plot.ts(v, main = "moving average")
Autoregesivo
w = rnorm(550,0,1)
x = filter (w, filter = c(1,-0.9), method = "recursive")[-(1:50)]
plot.ts(x, main = "autoregression")
Autoregesivo con drift
set.seed(154)
w = rnorm(200,0,1); x = cumsum (w)
wd = w + .2; xd = cumsum(wd)
plot.ts(xd, ylim = c(-5,55), main = "random walk")
lines(x); lines(0.2*(1:200), lty = "dashed")
ACF
par(mfrow = c(3,1))
acf(soi, 48, main = "Southern Oscilation Index")
acf(speech, 250)
par(mfrow=c(3,1))
acf(soi, 48, main="Southern Oscillation Index")
acf(rec, 48, main="Recruitment")
ccf(soi, rec, 48, main="SOI vs Recruitment", ylab="CCF")
# Varios vectores
persp(1:64, 1:36, soiltemp, phi=30, theta=30, scale=FALSE, expand=4, ticktype="detailed", xlab="rows", ylab="cols", zlab="temperature")
plot.ts(rowMeans(soiltemp), xlab="row", ylab="Average Temperature")
VECTORES AUTOREGRESIVOS
B1<-matrix(c(0.7, 0.2, 0.2, 0.7), 2)
var1<-VAR.sim(B=B1,n=100,include="none")
ts.plot(var1, type="l", col=c(1,2))
VARselect(var1, lag.max=9, type="const")## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 1 1 1 1
##
## $criteria
## 1 2 3 4 5 6
## AIC(n) -0.004781104 0.05373217 0.05239784 0.05727888 0.1294645 0.1453231
## HQ(n) 0.062008485 0.16504815 0.20824021 0.25764765 0.3743596 0.4347447
## SC(n) 0.160770072 0.32965079 0.43868392 0.55393241 0.7364854 0.8627116
## FPE(n) 0.995277884 1.05543577 1.05443691 1.06032595 1.1409261 1.1609733
## 7 8 9
## AIC(n) 0.1819582 0.2581421 0.3365710
## HQ(n) 0.5159061 0.6366164 0.7595717
## SC(n) 1.0097141 1.1962654 1.3850617
## FPE(n) 1.2068686 1.3060686 1.4176983var1_est <- VAR(var1,p = 1)## Warning in VAR(var1, p = 1): No column names supplied in y, using: y1, y2 , instead.var1_est##
## VAR Estimation Results:
## =======================
##
## Estimated coefficients for equation y1:
## =======================================
## Call:
## y1 = y1.l1 + y2.l1 + const
##
## y1.l1 y2.l1 const
## 0.67774002 0.23472999 -0.09862114
##
##
## Estimated coefficients for equation y2:
## =======================================
## Call:
## y2 = y1.l1 + y2.l1 + const
##
## y1.l1 y2.l1 const
## 0.41445790 0.58815891 0.01955353MODELOS ESTADO ESPACIO
El paquete base de R permite evaluar modelos estado espacio con la función link StructTS
# rm(list = ls())
fitNile <- StructTS(Nile, "level")
fitNile##
## Call:
## StructTS(x = Nile, type = "level")
##
## Variances:
## level epsilon
## 1469 15099tsdiag(fitNile)
plot(Nile, type = "l")
lines(fitted(fitNile), lty = "dashed", lwd = 2)
lines(tsSmooth(fitNile), lty = "dotted", lwd = 2)
SP_fit <- StructTS(Nile, "level")
SP_forecast <- predict(SP_fit, n.ahead = 100)$pred
SP_forecast_se <- predict(SP_fit, n.ahead = 100)$se
plot(Nile, type = "l")
points(SP_forecast, type = "l", col = 2)
points(SP_forecast - 2*SP_forecast_se, type = "l", col = 2, lty = 2)
points(SP_forecast + 2*SP_forecast_se, type = "l", col = 2, lty = 2)